/* Programa para testar a expressão da transformada de Fourier da biblioteca gsl, e a unidades corretas desta transformação */

/* Gerar uma função Gaussiana, e calcular a transformada de Fourier */

/* 	Calcular a transformada de Fourier da derivada da Gaussiana, com dois métodos: Calculando analiticamente a função (derivada da
Gaussiana) e depois calculando a transformada de Fourier, e depois calcular a transformada usando a expressão formal da transformada da
derivada, qual seja, iq vezes a transformada da função original. As duas transformadas devem prover o mesmo valor. */

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <complex.h>
#include <string.h>
#include <unistd.h>

#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_vector_complex.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_fft_complex.h>


/* gauss = \exp\left(-\frac{x^{2}}{2\gamma^{2}}\right) */
/* der_gauss = -\frac{x\exp\left(-\frac{x^{2}}{2\gamma^{2}}\right)}{\gamma^{2}} */


int main(void)
{
	int i, j;
	int N;
	int zetamin, zetamax, delzeta, kmin;
	double Dzeta, Dk;
	double gama;
	double tmp1, tmp2;
	double k;
	
	double *zeta;
	double *gauss;
	double *der_gauss;
	double *gaussK;
	double *der_gaussK;
	
	/* Os parâmetros N e delzeta são os importantes. Múltiplos de 2 somente, para facilitar */
	N = 1024;
	delzeta = 8;
	
	zetamin = -N/delzeta;
	zetamax = N/delzeta;
	// Dzeta = (zetamax - zetamin)/N;
	Dzeta = 2.0/delzeta;
	
	// kmin = -1/2*Dzeta 
	kmin = -delzeta/4;
	Dk = 1.0*delzeta/(2.0*N);
	
	gama = 10.0;
	
	zeta = malloc(N*sizeof(double));
	gauss = malloc(2*N*sizeof(double));
	gaussK = malloc(2*N*sizeof(double));
	der_gauss = malloc(2*N*sizeof(double));
	der_gaussK = malloc(2*N*sizeof(double));
		
	for ( i = 0 ; i < N ; i++ )
	{
		zeta[i] = zetamin + i*Dzeta;
		gauss[2*i] = exp((-1.0*zeta[i]*zeta[i])/(2.0*gama*gama)); /* Parte Real */
		gauss[2*i+1] = 0.0; /* Parte Imaginária */
		der_gauss[2*i] = -(zeta[i]/(gama*gama))*exp((-1.0*zeta[i]*zeta[i])/(2.0*gama*gama)); /* Parte Real */
		der_gauss[2*i+1] = 0.0; /* Parte Imaginária */
	}
	
	memcpy(gaussK, gauss, 2*N*sizeof(double));
	memcpy(der_gaussK, der_gauss, 2*N*sizeof(double));
	
	gsl_fft_complex_radix2_forward(gaussK, 1, N);
	gsl_fft_complex_radix2_forward(der_gaussK, 1, N);
	
	for ( i = 0 ; i < N/2 ; i++ )
	{
		j = i + N/2;
		k = (kmin + j*Dk);
		tmp1 = -2.0*M_PI*k*gaussK[2*i+1];
		tmp2 = 2.0*M_PI*k*gaussK[2*i];
		gaussK[2*i] = tmp1;
		gaussK[2*i+1] = tmp2;
	}
	for ( i = N/2 ; i < N ; i++ )
	{
		j = i - N/2;
		k = (kmin + j*Dk);
		tmp1 = -2.0*M_PI*k*gaussK[2*i+1];
		tmp2 = 2.0*M_PI*k*gaussK[2*i];
		gaussK[2*i] = tmp1;
		gaussK[2*i+1] = tmp2;
	}
	
	// imprimir
	for ( i = N/2 ; i < N ; i++ )
	{
		j = i - N/2;
		k = (kmin + j*Dk);
		printf ("%14.6f%14.6f%14.6f\n", k, gaussK[2*i]*gaussK[2*i] + gaussK[2*i+1]*gaussK[2*i+1], der_gaussK[2*i]*der_gaussK[2*i] +
		der_gaussK[2*i+1]*der_gaussK[2*i+1] );
	}
	for ( i = 0 ; i < N/2 ; i++ )
	{
		j = i + N/2;
		k = (kmin + j*Dk);
		printf ("%14.6f%14.6f%14.6f\n", k, gaussK[2*i]*gaussK[2*i] + gaussK[2*i+1]*gaussK[2*i+1], der_gaussK[2*i]*der_gaussK[2*i] +
		der_gaussK[2*i+1]*der_gaussK[2*i+1] );
	}

	return 0;
}

// gcc -Wall -lm -O3 -lgsl -lgslcblas teste_trans_fourier_gsl.c -o teste_trans_fourier_gsl.exe
